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ABSTRACT 

We present temporal and spectral analysis of simultaneous 0.5-79 keV Suiift-XKT and NuSTAR 
observations of the magnetar 4U 0142+61. The pulse profile changes significantly with photon energy 
between 3 and 35keV. The pulse fraction increases with energy, reaching a value of «20%, similar to 
that observed in IE 1841—045 and much lower than the «80% pulse fraction observed in IE 2259+586. 
We do not detect the 55-ks phase modulation reported in previous Suzaku -HXD observations. The 
phase-averaged spectrum of 4U 0142+61 above 20 keV is dominated by a hard power law with a photon 
index T H ~ 0.65, and the spectrum below 20keV can be described by two blackbodies, a blackbody 
plus a soft power law, or by a Comptonized blackbody model. We study the full phase-resolved spectra 
using the e ± outflow model of Beloborodov (2013). O ur results are c onsist ent with the parameters of 
the active j-bundle derived from INTEGRAL data by Hascoet et al. (20141. We find that a significant 
degeneracy appears in the inferred parameters if the footprint of the j-bundle is allowed to be a thin 
ring instead of a polar cap. The degeneracy is reduced when the footprint is required to be the hot 
spot inferred from the soft X-ray data. 

Subject headings: pulsars: individual (4LT 0142+61) - stars: magnetars - stars: neutron 


1. INTRODUCTION 


Magnetars — isolated neutron stars with inferred sur¬ 
face dipolar magnetic field strength, B SUI { > 10 14 G - 
were proposed to explain Soft Gamma Repeaters (SGRs) 
and later extended to include Ano malou s X-ra y Pul¬ 
sars (AXPs) (Thompson & Duncan 1995 1996)). Un¬ 
like canonical pulsars powered by their rotational en¬ 
ergy, the dominant energy reservoir of magnetars is their 
magnetic energy. The magnetar model attributes the 
anomalously high X-ray luminosity to the heating of 
the crust and magnetosphere by the dissipation of mag¬ 
netic energy. There are 28 magnetars that have been 
discovered to date, including 5 candidates sug gested o n 


the basis of the detection of X-ray bursts (Olausen & 
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KaspibOLflF 3 N otably, two magnetars, SG R 0418+5729 
QRea et al.||2(J10 ) and Swift J1822.3—1606 (Scholz et al. 


2014 and references therein) were recently shown to have 
-B S urf ~ 10 12 ~ 13 G, making them exceptions to the canon¬ 
ical classification above. For recent r eviews of ma gnetar 
obs ervations we refer the r eaders to Mereghetti (20081 
and Rea & Esposito (2011). 


The X-ray to y-ray spectrum of magnetars shows two 
peaks: a low-energy peak at ~0.5keV and a hig h-energy 
peak at energies greater than ^100 keV (see Kuiper et al. 
2006 Enoto et al.|2010 and references therein). The soft 
X-ray peak resembles a simple blackbody, likely originat¬ 
ing from the magnetar surface, with a tail extending to 
~10keV caused by radiative transfer of photons through 
the magnet ospheric plasma (see for example Thompson 
et al. 2002). Above energies of ~10-20kev, the hard 
X-ray component is dominant. 

The hard X-ray emission must be produced by the 
magnetosphere of the neutron star. A model making 
specific predictions for phase-resolved hard X -ray spec¬ 
tra emerged recently (Beloborodov] |2013b|a|) and was 
used to successfully fit the observations of IE 1841— 045, 
4U 0142+61, 1RXS J1708—4009, and IE 2259+586 (|Ah 


et al.||2013| |2015| |Hascoet et al.||2014| |Vogel et al.|20 

T he Nuclear Spectroscopic Telescope Array (NuST 


4 ). 


'The Nuc lear Spectroscopic Telescope Array ( NuSlAR ) 

(Harrison et al. 2013) with its excellent spectral and 
timing capabilities is well-suited for phase-resol ved spec¬ 
trosco py of magnetars in the 3-79 keV band (see) An et al. 
2014b for a review). In addition to 4U0142+61 AitW 
TAR has observed SGR 1745— 2900_J|MorietaL| 2013 


Kaspi et al. 2014 1, IE 1841—045 (An et al. 2013 ~|2015|) 


IE 1048.1-593 7 (lAn et al.|2014aD, and IE 2259+58BlfVc7 
gel et al.||2014 ). 


13 See the McGill Online Magnetar Catalog 
for a compilation of known magnetar properties: 
http: / / www.physics.mcgill.ca/^pulsar / magnetar/main, ht ml 
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1.1. 4U0142+61 

4U 0142+61 was discovered as a soft spectrum X-ray 
source in the Uhuru all sky survey, initially reporte d in 


the a nalysis of the first 70 days of data (Giacconi et al. 


1972). It remained an unexceptional source until 8.7-s 


X-rav pulsa tions were discovered by ASCA (Israel et al. 


1993 1994). 

The soft X-ray emission from 4U 0142+61 is well de¬ 


scribed bv a blackbodv with k^T ~ 0.4 keV and a power 

law with index T ~ 3.7 (White et al. 1996 

Israel et al. 

1999; Paul et al. 2000; Juett et al. 2002; Pate 

et al. 2003; 

Gohier et al. 2004, 2005 Rea et al. 2007 

3noto et al. 

2011). Pulsed high-energy emission was c 

letected be- 


instrument on INTEGRAL (den Hartog et al. 2004). The 
hard X-ray spectrum (>10kcV) is domina ted by a power- 


law component with F ~ 1 dKuiper et al. 2006 

den Har- 

tog et al. 

2006 

2008 bj Enoto et al. 

T01T). Using up- 


telescopes, the hard X-ray power-law cutoff energy was 
sugge sted to be between ~ 200-750 keV pvuiper et al. 
2006). 


The soft X-ray pulse profiles of 4U 0142+61 were shown 
to undergo long-term changes using RXTE observations 
spread over 10 years (Dib et al. 2007) and l ater Chandra , 
XMM -Newton and Swift observations (Gonzalez et al. 
20101. The pulse fractions were observed to increase over 
time leading up to a grou p of three burst s that occured 
between 2006-2007 ( Gonzalez et al.pOl O). 

There has been much debate about the intrinsic soft 
X-ray spectra of magnetars, the measurement of which 
depends o n the absorptio n column, Ah, along the line 
of sight. Durant & van Kerkwijk (2006b), hereafter 
D06b, estimated the A H to 4U 0142+61 to be (6.4 ± 
0.7) x 10 21 cm -2 by fitting high-resolution grating spec¬ 
tra around individual photoelectric absorption edges of 
oxygen, iron, neon, magnesium and silicon. They also 
showed that the abundance ratios of Ne/Mg and O/Mg 
fo r 4U 0142+61 ar e closer to the revised solar abundances 
of lAsplund et al.| fl2005|) compared t o the old standard 
abundances of Anders & Grevesse (1989). This mea¬ 
surement has the advantage of being less sensitive to the 
choice of model used to describe the intrinsic magne- 
tar spectrum. However, since only the data near photo¬ 
electric edges is fitted, the fitting requires high-quality 
X-ray data. 

Unlike the IVh measurements from the high-resolution 
X-ray spectra, most fits of the low-energy spectrum (« 
0.5-10 keV) with a blackbody pl us power-l aw model con¬ 
verge to IVh ~ 1.0 x 10 22 cm~ 2 ( Rea et aL||2007 and ref¬ 
erences therein) w hich used the a bun dances from Anders 
& Grevesse| (|1989D . We note that |Enoto et al.| (|2UIl|) ob¬ 
tained IVh values consistent with the D06b value from 
broadband spectral fits to Suzaku data, however, no 
abundance model was specifi ed. In this work, we use 
solar abundance values from Asplund et al. (2009) - 
the ‘aspl’ model in XSPEC — as default, and we also test 
our fits with other abundance models. 

By identifying core helium-burning giant stars - i.e. 
red clump stars — from the 2MASS catalog and es¬ 
timating the variation of optical extinction as a func¬ 
tion +>f distance in t he direction of magnetars, |Du-| 
rant & van Kerkwijk (2006a) estimated the distance to 


4U 0142+61 to be 3.6 ± 0.4 kpc. This distance estimate 
used the Ah = (6.4 ± 0.7) x 10 21 cm -2 estimated from 
the photo-electric absorption edges. Previous measure¬ 
ments of optical extinction (Ay) have concluded that 


the Ay for 4U 0142+61 should be less than 5 (Hulleman 


et al. 2004 

Durant & van Kerkwijk 2006bl), correspond- 

ing to TVh < 9 x 10 2i cm 2 ( 

Predehl & Schmitt 

1995 


Ay = TVh x 5.6 x 10“ cm"). 

In this paper, we present a phase-resolved spectral and 
timing analysis of coordinated Swi ft-XRT and NuSTAR 
spectra of 4U0142+61. We use the|Hascoet et al.|([2bl41 
framework to test the electron-positron outflow model 
and constrain physical parameters. The paper is orga¬ 
nized as follows. In Section [2j we describe our observa¬ 
tions, data and data reduction procedure. In Section [3] 
we describe the results of our timing analysis, spectral 
analysis and model fitting. In Section |4j we discuss 
these results in the context of previous observations of 
4U 0142+61 and those of other magnetars. 

2. OBSERVATIONS AND ANALYSIS 


NuSTAR (Harrison et al. 2013) is a 3-79keV focus¬ 
ing hard X-ray mission. It consists of two identical co¬ 
aligned Wolter-I telescopes with CdZnTe detectors at the 
focal planes. The telescopes provide a point-spread func¬ 
tion with a half-power diameter (HPD) of 58" over a 
field of view of 12'xl2'. The energy resolution varies 
from 0.4 keV at 6keV to 0.9 keV at 60keV. The two fo¬ 
cal plane modules are referred to as FPMA and FPMB. 
4U 0142+61 was observed by NuSTAR between 2014 Mar 
27 30 during a 44-ks observation simultaneous with a 

24-k s observatio n with the Swift X-ray Telescope (XRT; 
Burrows et al.|2005|). The details of the observations are 
summarized in Table [Tj 

We performed the processing and filtering of the NuS¬ 
TAR event data with the standard NuSTAR pipeline 
version 1.4.1 and HEAS0FT version 6.16. We used the 
barycorr tool to correct the photon arrival times for the 
orbital motion of the satellite and the Earth at the opti¬ 
cal position of 4U0142+61 — a = 01 h 4 6 m 22 s .407, £ = 
+61°4 5'03".19 (J2000) — as reported by Hulleman et al. 
(2004). The source events were extracted within a 50- 
pixel (120") radius around the centroid and suitable 
background regions were used. Spectra were extracted 
using the mrproducts script. Using grppha, all pho¬ 
tons below channel 35 (3keV) and above channel 1935 
(79 keV) were flagged as bad and all good photons were 
binned in energy to achieve a minimum of 30 photons 
per bin. 

The Swift-XRT data were obtained in the Windowed 
Timing (WT) mode and were processed with the standard 
xrtpipeline and the photon arrival times were corrected 
using barycorr. The xrtproducts script was used to ex¬ 
tract spectra and lightcurves within a radius of 25 pixels 
(59")[2 Photons in channels 0-29 (energy < 0.3 keV) 
were ignored and all channels between 0.3-10 keV were 
binned to ensure a minimum of 30 photons per bin. 

The Swift-XRT and NuSTAR FPMA and FP MB spec¬ 
tra w ere fit simultaneously in XSPEC vl2.8.1 (Arnaud 


1996) using two freely varying cross-normalization con¬ 
stants, assuming the normalization of Swift-XRT to 

14 As per the Swijl.-XUT data analysis thread: http://www. 
swift.ac.uk/analysis/xrt/spectra.php 
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Table 1 

X-ray observations of 4U 0142+61 in 2014. 


Obs ID 

Start 

(UT) 

End 

(UT) 

Exp 

(ks) 

Rate a 

cts/s 

NuSTAR 





30001023002 

Mar 27 13:35 

Mar 28 00:45 

7 

1.3 

30001023003 

Mar 28 00:45 

Mar 30 13:00 

37 

1.3 

Swift-X RT (Windowed Timing Mode) 



00080026001 

Mar 27 13:36 

Mar 27 21:52 

4.9 

4.2 

00080026002 

Mar 28 07:10 

Mar 29 23:15 

12.9 

4.0 

00080026003 

Mar 30 00:42 

Mar 30 08:57 

6.6 

4.1 

a 0.5-10 keV count rate for Swift-XRT and 3- 

-79 keV count 


rate from NuSTAR. 


servations of 


Gonzalez et al. 


(2010) from 2008 March. 


However, the 20-35 ke V and 35-50 keV observations from 
NuSTAR show a double peak structure with the primary 
peak having approximately twice the peak amplitude as 
compared to the secondary peak. The corresponding IN¬ 
TEGRAL pulse profile reported in dH08 showed a double 
peak structure with both peaks of equal amplitude. This 
difference is also present in the pulse fraction analysis 


presented in Section 3.3 


3.2. Pulse Morphology 

To explore the variation in pulse shape as a function of 
energy, we decomposed the pulses into Fourier harmon¬ 
ics. We define Fourier coefficients ak and bk as 


be fixed to unity. Timing analysis was performed on 
exposure-corrected lightcurves and event lists using cus¬ 
tom MATLAB scripts. 



N 

X^ cos 

1=1 



and 


( 1 ) 


3. RESULTS 
3.1. Pulse Profile 

We analyzed the barycentered 3-79 keV NuSTAR 
events using epoch folding (|Leahy| 1987) and mea¬ 
sured the rotation period of 4U 0142+61 to be P = 
8.689158(4) s. This is consistent with the period mea¬ 
sured with the Swift-XRT observation and is also consis¬ 
tent with the period P = 8.689163(5) s expected at the 
epoch of observation based on the last ephemeris mea¬ 
sured after the glitch of 2011 July, reported by IDib &| 
Kaspi (20l4|. 

We folded the Swift-XRT and NuSTAR events in 8 
energy bands into 20 phase bins with our measured pe¬ 
riod to compare the pulse morphology as a function of 
energy (Figure!!]). The energy bands 0.3-1.5keV, 1.5- 
3 keV, 3-5 keV, 5-8 keV, 8-20 keV, 20-35 keV, 35-50 keV 
and 50-79 keV — were chosen to have approximately 
equal counts in each band. The 3-5 keV and 5-8 keV 
data from Swift-XRT had far lower count rates than the 
corresponding NuSTAR observations and hence were not 
used for the final analysis. However, we confirmed that 
the NuSTAR and Swift-XRT pulse profiles in these two 
overlapping energy bands are consistent within the error 
bars. 

There is a clear gradual change in the pulse morphol¬ 
ogy as the energy band crosses ^3keV and ~20keV, cor¬ 
responding to the different spectral components — mod¬ 
ified blackbody or hard power law — that dominate the 
spectrum at these energies. This change in morphology is 
also observed in the dom inance of the Fourier harmonics 
described in Section [+2l The 0.3-1.5keV and 1.5-3 keV 
pulse profiles consist of two peaks at phases of f = 0.3 
and 4> = 0.6 separated by a sharp dip at f> = 0.5. In 
the 3-5 keV and 5-8 keV bands, the peak at <j> = 0.3 
dominates the pulse and the dip at <f = 0.5 deepens sig¬ 
nificantly. There is a small dip at <f> = 0.9 separating a 
possible second pulse peak from the primary. Moving to 
higher photon energies, in the 8-20 keV and 20-35 keV 
bands, a second pulse rises in amplitude at <fi ~ 0.65 to¬ 
wards energies of 50 keV. The primary pulse also shows 
signs of broadening as a function of energy. 

Our low-energy result s are consistent w i th the RXTE 
observations reported by den Hartog et al. (2008b) (here¬ 
after dH08) and with the 0.5-10 keV XMM-Newton ob- 



where N is the number of phase bins and pj is the num¬ 
ber of photons in each phase bin and j and k are indices 
referring to the phase bins and the Fourier harmonics re¬ 
spectively. We define the strength of each Fourier com¬ 
ponent to be Ak = \/a\ + b\. We define A to ta\ as 


dtotal 


\ 


N 




( 3 ) 


k= 1 


We find that most of the variational power in the pulses 
is explained in the first six harmonic coefficients. The dis¬ 
tinct variation of pulse shapes with energy can be seen 
in Figure [2] The fraction of power in the first harmonic 
(Ai/H to t a i) decreases with energy until approximately 
40 keV and then increases, whereas the fraction of power 
in the second harmonic (A 2 /A tota i) increases with energy 
until approximately 40 keV and then decreases. 

This behavior of the harmonics is significantly dif¬ 
ferent from that of IE 2259+586 presented in Vogel| 
et al. (2014). In IE2259+586, the normalized Ai value 


increases as a function of energy until approximately 
12keV and then decreases as a function of energy. The 
normalized value of H 2 decreases as a function of energy 
until approximately 12-15 keV and then increases. 


3.3. Pulse Fraction 

We quantify the strength of the pulsations using two 
different methods. We define the root-mean-square pulse 
fraction as 


PFrms 


+ -«+<)) 

a 0 


( 4 ) 


where ak and bk are as defined above and cr ak and <Jb k 
are the uncertainties in ak and bk , respectively, calculated 
using Poisson variances as 


a 


2 

ak 


1 

7V2 


E4 


3 = 1 


cos 2 



( 5 ) 
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Rotation Phase 



Rotation Phase 



Rotation Phase 



Figure 1. Swift-XKT and NuSTAR pulse profiles in different energy bands. The annotation in the upper left corner of each plot specifies 
the telescope (‘S’: Swift-X RT, ‘N’: NuSTAR) and the energy band for each plot. The last plot is the total 0.3-79keV (Swift-X RT and 
NuSTAR) count rate ( norm alized to the average value) marked with phase bins — ‘A’, ‘B’ and ‘C’ and ‘DC’ — used for fitting the e~ — e~*~ 
outflow model (Section |3.7| ). Two pulse periods are shown for clarity. 

N 




= Iy 

N 2 


N 2 


9 • 2 

ai .. sin 

P 3 


3 =1 




( 6 ) 


This definition, including the correction term, a 2 k + a 2 k , 
has been shown to be a robust and accurate metric of 


p ulse fraction in noisy data (see Appendix 1 of An et al. 


2015 for a detailed discussion). 


We also define the area pulse fraction described by 


Gonzalez et al. (2010) as 


EyLi Pj 

\~~\N 

£,=i Pi 


PF — 

- L -L a.rfia. — 


(7) 


This definition is consistent with that used by dH08. 
However, it is challenging to determine the true value 
of min(pj), and both noise and binning tend t o bias the 
PTkrea metric upwards by as much as 20% (An et al. 


20151. 

figure [3] shows the variation of PF aiea (filled symbols) 
and PFrms (empty symbols) as a function of energy. 
Note that while our measurements of PF area have an 
increasing trend at energies >10keV, the PF area values 


above 20 keV are also consistent with a constant value of 
«35%. The near-linear increase in PP area as a function 
of energy is consistent with the results of dH08, though 
we note that our PF area measurements are consistent ly 
higher than those of dH08 and those of |Gonzalez et ah] 
(20101. The RMS pulse fraction, PFrmS) increases with 
energy up to an energy of 35 keV. However, the absolute 
normalization is different due to the different definitions 
of pulse fractions. The possible decrease in PFrms in 
the 35-50 keV and 50-79 keV bands may be due to the 
emergence of two nearly equal amplitude peaks in the 
pulse profile with lower count rates. A similar reduction 
in RMS pulsed fraction with energy had been reported 


for IE 1841—045 in the 16-24 keV band observations (An 
et al.|2013). However, with more NuSTAR observations 


the variations were s hown to be dep endent on the ex¬ 
act energy bins used (An et al.||20151. The overall trend 
of both RMS and area pulse fraction was shown to in¬ 
crease with energy, with PFrms increasing up to 20% at 
50 keV and PF aiea increasing to 50% at 50 keV. Similar 
to IE 1841—045, PFrms does not show signs of increas¬ 
ing to 100% with increasing energy as was suggested from 
the INTEGRAL data (dH08). 
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Figure 2. Variation in the first (filled symbols) and second (empty 
symbols) harmonic amplitudes as a function of photon energy. 
Both values are normalized with respect to the total amplitude 
of the variation (A tota i). The Swift-X RT data points are shown as 
red circles and the NuSTAR data points are shown as black squares. 
The filled areas (solid for Ai/A tota i and hashed for A2/A to tai) show 
the l-cr error regions. (A color version of this figure is available in 
the online journal.) 



Figure 3. Variation in the area pulsed fraction (filled symbols) 
and root-mean-square pulsed fraction (empty symbols) as a func¬ 
tion of energ y. The Swift-X RT and NuSTAR symbols are the same 
as in Figure [2] The filled areas (solid for PF area and hashed for 
PFjims) show the l-cr error regions. A color version of this figure 
is available in the online journal. 


3.4. Non-Detection of Precession 


Makishima et al. (2014) (hereafter ME14) reported 
a phase modulation m the 8.7-s rotation period of 
4U 0142+61 with an amplitude of 0.7 s and a period 
of 55±4ks (ssl5hr) detected from 15-40 keV HXD-PIN 
data gathered with Suzaku in 2009 August. This was 
interpreted as possible evidence for the precession of the 
neutron star caused by slight deviation from spherical 
symmetry. The same search in 5it^aA:M-HXD data gath¬ 
ered in 2007 August and XIS data from 2007 August 
and 2009 August did not lead to detection of preces¬ 
sion. Since our observations are spread over 4 days, we 
searched for possible variations in the rotation perio d or 
rotation phase following the same Z „ analysis (Brazier 



- 0.0020 - 0.0015 - 0.0010 - 0.0005 0.0000 0.0005 0.0010 0.0015 0.002 

Period (s) - 8.689158 


Figure 4. Variation of +( as a function of rotation period for 
4U 0142+61 for 15-40 keV band using NuSTAR data (solid black 
curve) overlaid on Suzaku -HXD results from ME14 (red dashed 
and solid curves). We find that Z\ peaks to a va lue o f Ri360 at the 
measured rotation period of 8.689158 s (Section . The dotted 
red line is Z\ for the raw (non-demodulated) Suzaku -HXD data 
reported in Figure lb of ME14. The solid red line is for the 
same HXD data after optimal demodulation. The vertical dashed 
red line is the high-energy (15-40 keV) rotational period reported 
by ME14 with the shaded region denoting the reported error and 
the vertical dashed red line the low-energy (0.3-10 keV) rotational 
period reported from the XIS data for the same epoch. Vertical 
black lines are the expected rotation periods for the 2009 August 
and 2014 March epochs from the RXTE data. The rotation periods 
marks are offset vertically for clarity. 


19941 steps reported by MEM. 


For n = 3 and n = 4, we find that Z‘* peaks to a value 
of «360 at the rotation period of P = 8.689158(4) s for 
the NuSTAR data without the need for demodulation 
(FigureEj). This is consistent with the P = 8.689163( 5) s 
expected * from the RXTE ephemeris of 4U0142+61(Dib 


& Kaspil 2014 Ephemeris E, Table 6). From the 2009 


August 12-14, Suzaku-XIS data, ME14 reported a period 
of P = 8.68891 +0.00010 s, which is inconsistent with the 
value of P = 8.68869734(8) s r eported from the RXTE 
ephemeris (|Dib & Kaspi; 2014 Ephemeris D, Table 6). 
The Suzaku-XIS and HXD measured rotation periods are 
marked in Figure [4] with their corresponding errors along 
with the RXTE ephemeris for comparison. 

Figure [4] indicates that in the NuSTAR observations, 
the pulsations were detected at a higher significance than 
during the high-energy Suzaku observations. This result 
is (a) similar in value and shape to the result reported 
in Figure la of ME14 (XIS data from the same Suzaku 
observations), (b) significantly higher than the Z\ = 12 
and Z\ = 16 (without demodulation) and Z\ ss 52 (after 
optimal demodulation) reported in their Figures lb and 
lc. 


We searched for phase modulation in the data by 
shifting the arrival times of each photon by At = 
Asm(2nt/T — (f> o), where t is the time of arrival, A is 
the modulation amplitude (with units of time), T is the 
modulation period and <f>o is the initial phase. We mea¬ 
sured Z\ after varying T between 45 ks to 65 ks in steps of 
2.5 ks, A between 0-1.2 s in steps of 0.1 s and (/>o between 
0°-360° in steps of 20°. These results were compared to 
Figure 2 of ME14. We find that unlike ME14, Z\ peaks 
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to a value >350 at A = 0 reducing to ~100 at A = 1.2. 
Z\ is also nearly independent of <j >o at any given T and 
A. We find no preference for the values of A = 0.7 + 0.3s 
and 0 O = 75° ± 30°. 

Assuming a 55-ks period reported by ME14, we split 
the data into six subsets, each 9.17 ks long, and cre¬ 
ated individual pulse profiles by folding each subset 
at P = 8.689158(4) s. We find no phase change be¬ 
tween any two pulse profiles (compared to Figure 3 of 
ME14). We also find that the post-demodulation pulse 
profile reported in Figure If of ME14 is triple-peaked 
with PF aiea < 10% and significantly different from the 
double-peaked NuSTAR profile and the 20-50 keV pulse 
profile reported from XMM-Newton data in dH08 (Fig¬ 
ure 7E), each with PF area « 30%. 


3.5. Spectral Fits 

The phase-averaged X-ray spectrum of 4U 0142+61 has 
been previously fit with a hard power law at high ener- 
gies (>20keV) and by a modified blackblody (BB) or 
combination of blackbodies at low energies (<10keV). 
We fit the extracted spectrum in XSPEC with three dif¬ 
ferent models: I) a hard high-energy power law (PL) plus 
a blackbody and a soft low-energy PL, II) a hard high- 
energy PL plus two BB models at low energies, and III) 
a hard high -ener gy PL plus a comptonized BB (nthcomp 
Zycki et al. 1999) at l ow energies. Each model included 
a tbabs model (| Wilms et al.||2000|) with solar elemental 
abundances ( 1 aspl 1 ) and cross- sections described by the 


bcmc’ model (Balucinska-Church & McCammon 1992 


Yan et al. 19981) to fit for photo-electric absorption anc 


a cross-normalization parameter to allow for slight cal¬ 
ibration differences between the Swift-XRT , NuS TAR 
FPMA and NuSTAR FPMB detectors. In Section |3.7[ 
we present fits to the phase-averaged spectra with a cus¬ 
tomized combination of model I and model II: a sum 
of one blackbody, one modified blackbody with a soft 
power-law tail and a hard power law, si milar t o th e reso - 
nant compton scattering model used by Rea et al. (2007). 

The results of the fitting are shown in Table [2] We 
checked the validity of each model fit by generating 1000 
sets of synthetic data based on the best-fit model param¬ 
eters and testing their % 2 with respect to the model. If 
the distribution of y 2 values from synthetic data is sig¬ 
nificantly lower than the x 2 of the real data, the fit is 
deemed to be unacceptable. The ‘goodness’ parameter 
in Table |2] shows the fraction of synthetic y 2 values that 
are lower than the \ 2 from the real data. Model I and II 
have traditionally been used to describe the spectrum of 
4U 0142+61. However, we find that while the models can 
match the spectral distribution visually, the fits are sta¬ 
tistically unacceptable. Model III provides a statistically 
acceptable fit. 

Figure [ 5 ] shows the fit of model I (BB+2PL, top 
panel), model II (2BB+PL, middle panel) and model III 
(nthcomp+PL, bottom panel). As noted in Table [2j the 
index of the high-energy power law (T h) varies between 
0.3-1.0 depending on the model used to fit the low en¬ 
ergy spectrum. This is reflected in the residuals at the 
high-energy end of Figure [5j It is clear that the best- 
fit high-energy power law underpredicts the data at high 
energy for model II (middle panel) while it over predicts 
the data when the < lOkeV spectrum is modeled with 


model I (top panel). Note that while the nthcomp model 
is a good phenomenological fit to the low and intermedi¬ 
ate energy X-ray spectrum and the blackbody emission 
from the surface is expected to be upscattered by high- 
energy electrons outside the neutron star, the nthcomp 
model does not accurately account for the effects of the 
extremely strong magnetic field on the photon scattering 
process. 

If we restrict the fits to the low-energy spectrum (< 
lOkeV), model I fits improve with x 2 /dof = 1554.3/1391 
and the parameter values are similar (within 2-cr) to 
those in Tabic [2j suggesting that this model can fit the 
low-energy spectrum well but cannot describe the 10— 
20keV region of the spectrum. Fitting model II to the 
low-energy spectrum produces x 2 /dof = 2082.4/1391, 
which is statistically unacceptable. 

Assuming a nominal value for the neutron star radi us 


-Rns = 10km a nd a distance of 3.6kpc (Durant & van 


Kerkwijk 2006a), we can calculate the fraction of the neu¬ 


tron star surface area (.4ns) covered by the blackbody of 
a given flux normalization. For model I, the blackbody 
has a bolometric luminosity of Lboi = 1-3 x 10 35 ergs -1 , 
covering 0.2 Ans and contributing 11% of the 0.5-79 keV 
X-ray luminosity (most of it in the 0.5-10keV band), 
while the soft power law contributes the remaining 84%. 
For model II, the low temperature blackbody (Lboi = 
2.5 x 10 35 ergs -1 ) covers 0.6 Ans and contributes 75% 
of the luminosity, while the high temperature blackbody 
has a luminosity of Lboi = 3.3 x 10 34 ergs -1 emanat¬ 
ing from a hotspot covering 0.004 Ans of the surface and 
contributing 10% of the X-ray luminosity. In model III, 
the nthcomp component (combined blackbody and comp¬ 
tonized power law) contribute 86% of the total X-ray 
flux, similar to the contributions of the two blackbodies 
of model II. 


3.5.1. High-Energy Power Law 

The differences in the hard power law between different 
models are caused by the inability of these phenomeno¬ 
logical models to accurately describe the spectrum be¬ 
tween approximately 10-20 keV. In order to minimize the 
figure-of-merit (% 2 in this case) for the fit, XSPEC forces 
variations in the hard power-law index and normaliza¬ 
tion. To better measure the slope of the hard power 
law, we restricted the energy range from 20-79 keV and 
fit the phase-averaged spectrum with a power law. We 
measure Tjj = 0.65 ± 0.09. The 2 parameter fit yielded 
a Xred = 477.12 for 492 degrees of freedom and a p- 
value of 0.68. This is independent of TVh and the model 
used to describe the soft X-ray (< lOkeV) spectrum. 
When the energy range is further constrained (i.e. in 
the ranges 25-79keV and 30-79keV), we get consistent 
measures of T# but with larger uncertainties. This value 
is lower than the Tjj = 0.93 + 0.06 measured by dH08 
and Tfj = 0.89 ± 0.10 measured by|Enoto et al.| (|2011|). 
However, note that the values measured by dH08 varied 
from 0.79 ± 0.10 to 1.21 ± 0.16 over different datasets. 
The nthcomp model provides the least structured residu¬ 
als and the value of T# is closest to the high-energy-only 
value. 

If the high-energy hard power law is frozen to the 
T h = 0.65 value and the corresponding normalization 
and the low-energy (0.5-10 keV) spectrum is described 




























NuSTAR and Swift-XRT Observations of 4U 0142+61 


7 


Table 2 

Phase-averaged spectral fits. 


Component Parameter Value 


Model I 


const*tbabs*(bbody+powerlaw+powerlaw) 


const 

Cfpma 

0.981 ± 0.015 


Cfpmb 

0.977 ± 0.015 

tbabs 

N b (10 22 cm -2 ) 

1.30 ±0.03 

bbody 

k B T BB (keV) 

0.462 ± 0.005 


norm 3 ( 10 -3 ) 

1.06 ±0.04 


F b b c 

0.11 

powerlaw 

r s 

3.85 ±0.04 


norm b 

0.18 ± 0.01 


F B L,S C 

0.84 

powerlaw 

r H 

0.29 ±0.05 


norm b ( 10 -5 ) 

2.3 ±0.4 


^PL,H C 

0.05 

Xr 2 ed/ d ° f 

1.174/2408 

p- value 


4.8 x 10 “ 9 

goodness d 


100 % 

Model II 



const*tbabs* 

(bbody+bbody+powerlaw) 

const 

Cfpma 

1.033 ±0.016 


Cfpmb 

1.029 ±0.016 

tbabs 

AT h (10 22 cm -2 ) 

0.52 ±0.01 

bbody 

k B T BBt i (keV) 

0.422 ±0.004 


norm 3 ( 10 -3 ) 

1.90 ±0.02 


+BB,l c 

0.75 

bbody 

k B 'T BB '2 (keV) 

0.93 ±0.02 


norm 3 ( 10 -4 ) 

2.5 ±0.1 


+BB,2 C 

0.10 

powerlaw 

V H 

1.03 ±0.05 


norm b ( 10 -4 ) 

9 7+0.4 
' —0.3 


Cpl,h c 

0.15 

Xr 2 ed/ d ° f 


1.130/2408 

p- value 


6.7 x 10 _b 

goodness d 


99.7% 

Model III 



const*tbabs* 

(nthcomp+powerlaw) 

const 

Cfpma 

1.001 ± 0.015 


Cfpmb 

0.998 ± 0.015 

tbabs 

ATh ( 10 22 cm -2 ) 

0.65 ±0.02 

nthcomp 

r s 

4.86 ±0.04 


k B T BB (keV) 

0.346 ± 0.004 


k B T e - (keV) 

> 37.3 


norm b ( 10 " 2 ) 

6.5 ±0.2 


fP c 

r nthcomp 

0.86 

powerlaw 

V H 

u - ‘°-0.04 


norm b ( 10 -4 ) 

1 1+02 
1 * 1 -0.1 


F bb c 

0.14 

Xr 2 ed/ d ° f 


1.07/2408 

p -value 


6.3 x 10 " 3 

goodness d 


82.8% 


a Normalization in units of L^q/D^q, where L 39 
is the source luminosity in units of 10 39 ergs _1 
and Dio is the distance to the source in units of 
10 kpc. 

b Normalization in units of 
photons keV —1 cm -2 s -1 at 1 keV. 
c Fraction of the total 0.5-79 keV flux contributed 
by the component. 

d Goodness of fit is the percentage of x 2 values 
from 1000 Monte Carlo simulations synthesized 
the best-fit model parameters that are less than 
the best-fit x 2 value. The data are indistinguish¬ 
able from the synthesized data if the goodness 
«50%. 


1 



CO 


1 = 0.01 

o 

CO 

o 10- 3 

o 



1 



1 


10 


Energy (keV) 



Energy (keV) 



Energy (keV) 


Figure 5. Unfolded phase-averaged Swift-XRT and NuSTAR 
spectrum and the ratio of the data to the model. The model fit 
shown is const*tbabs* (bbody+powerlaw+powerlaw) (Model I, top 
panel), const*tbabs*(bbody+bbody+powerlaw) (Model II, middle 
panel) and const*tbabs*(nthcomp+powerlaw) (Model III, bottom 
panel). The colors are as follows: black: NuSTAR FPMA Obs I, 
red: NuSTAR FPMA Obs II, blue: NuSTAR FPMB Obs I, green: 
NuSTAR FPMB Obs II, cyan, yellow, magenta: Swift-XRT Obs I, 
II and III. 
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2 const*tbabs*(bbody+powerlaw+powerlaw) 



2 const*tbabs*(bbody+bbody+powerlaw) 



2 const*tbabs*(nthcomp+powerlaw) 



Energy (keV) 


Figure 6. Data to model ratio for models fit with Nn and hard 
power-law param eters frozen to independently measured values 
(see Section|3.5.2|). From top to bottom, the plots represent model 
I, model II and model III. 


with parameters from model I, the fit worsens with 
X 2 = 3104.0 in 2412 degrees of freedom. The fit for 
model II also worsens with x' 2 = 3061.3 in 2412 degrees 
of freedom. Both models show structured wavy residuals 
between 5-11 keV suggesting that the models are failing 
to capture all the structure in the data. Model III fit pa¬ 
rameters do not change values within 1-cr errors when the 
hard power law is frozen, y 2 = 2602.8 for 2412 degrees 
of freedom and the goodness of fit is 85%. 


3.5.2. Photoelectric Absorption 

We note that the 1V H value from Model III (Table [2]) is 
consistent with the D06b value of (6.4±0.7) x 10 21 cm -2 . 
To check on the influence of various abundance models, 
we re-tested the model fits with the seven abundance 
data sets used in X-ray astronomy (Table[3j in chronolog¬ 
ical order). We find that the older abundance data sets 
(aneb, angr) consistently provide lower Nn values com¬ 
pared to the newer data sets (wilm.lodd, aspl) and that 
the nthcomp+PL model provides the better fit irrespec¬ 
tive of which abundance model is used. Note that these 
valu es are roughly sim ilar to the previous values reported 
by |Patel et at fl2003|) (7V H = (0.93 ± 0.02) x 10 22 cm~ 2 


using aneb abundances and Model I) and by Rea et al. 


(2007) (IVh = (0.926 ± 0.005) x 10 22 cm , using angr 
abundances and Model I). 


3.5.3. Freezing Nu and the High-Energy PL 

The D06b value of IVh and our measurement of the 
high-energy PL are independent of the complicated spec¬ 
tral shape at low and intermediate energies (< 20keV). 
By freezing the value of Nb = 6.4 x 10 21 cm -2 and freez¬ 
ing the high-energy power la w to t he slope and normal¬ 
ization measured in Section |3.5.1[ we can explore the 
low-energy spectral shape and investigate whether addi¬ 
tional spectral components are required to fully describe 
the low-energy distribution. For technical reason^] 

15 Since the TVh value affects only the Swift-XRT spectrum and 
the high-energy PL affects only the NuSTAR spectrum, allowing 
the cross-normalization factors to vary freely effectively allows the 


the cross-normalization factors between Swift-XRT and 
NuSTAR were frozen to unity. 

We find that Model I and Model II fits worsen signif¬ 
icantly with x 2 = 5896.6 and y 2 = 3480.0 for 2413 de¬ 
grees of freedom respectively with extremely wavy resid¬ 
uals (Figure [6j Panels I and II). Model III provides a 
fit parameters similar to that from Table [2] with a total 
X 2 = 2608.3 for 2413 degrees of freedom. 


3.6. Phase-Resolved Spectral Fits 

We created good-time-interval (gti) files using the 
measured period of 4U 0142+61 and extracted Swift- 
XRT and NuSTAR spectra in five equal phase bins: 
f> =0.0-0.2, 0.2-0.4, 0.4-0.6, 0.6-0.8 and 0.8-1.0. We fit 
each 0.5-79 keV spectrum with a BB+2PL and nthcomp 
+ PL models. The fit parameters are detailed in Ta¬ 
ble [4] We froze the values of Cfpma, Cfpmb and -/V H to 
those fit in the phase-averaged spectrum using the same 
spectral model (as in Table [2]). 

The y 2 ed for each individual phase is lower than that 
from the corresponding fits of the phase-averaged spec¬ 
tra. Figure [7] shows fit parameters for the nthcomp+PL 
model as a function of phase compared to the phase- 
averaged fit values. The spectral shape parameters T#, 
IzbTbb and Fg are statistically consistent within 3-cr 
with the values measured from the phase-averaged spec¬ 
tra. However, we detect a very significant increase in 
the hard power-law normalization in the 0.2-0.4 phase 
range which corresponds to the peak of the high-energy 
pulse profiles (20-35 keV, 35-50 keV, Figure [l|. Simi¬ 
larly, the normalization of the nthcomp component shows 
a sharp decrease in the 0.4-0.6 phase bin which corre¬ 
sponds to the dip at (j) = 0.5 in the 3-5 keV and 5-8 keV 
pulse profiles. For the BB+2PL model, we observe a 
similar trend with the hard power-law normalization sig¬ 
nificantly increasing in the 0.2-0.4 phase range and the 
soft power-law normalization (which contributes approx¬ 
imately 85% of the X-ray flux ) de creasing between the 
0.4-0.6 phase range. In Section [3+1 we describe the vari¬ 
ation in the high-energy spectra in greater detail with a 
physical emission model. 


3.7. e ± Outflow Model 

Next we_test the coronal outflow model proposed by 
(2013a). The model envisions an outflow of 
ectron-positron (e ± ) pairs created by elec- 


Beloborodov 
relativistic e 


trie discharge near the neutron star. The outflow moves 
along the magnetic field lines and gradually decelerates 
as it (resonantly) scatters the thermal X-rays. The out¬ 
flow fills the active “j-bundle” that carries the elec¬ 
tric currents o f twisted magnetosplreric field lines (Ben 


loborodov 2009). It radiates most of its kinetic energy 
in hard X-rays before the e ± pairs reach the top of the 
twisted magnetic loop and annihilate. 

The magnetic dipole moment of 4U 0142+61 is p ss 
1.3 x 10 32 Gem 3 ( calculated fro m the spin-down rate ; 
Dib & Kaspi 2014). Similar to |Hascoet et al. (2014), 
we assume a simple geometry where the j -bundle is ax- 
isymmetric around the magnetic dipole axis. However, 


high-energy PL normalization to vary, spoiling the high-energy 
fit. To prevent this effect, we must freeze the cross-normalization 
constants. The expected systematic cross-calibrat ion error be¬ 
tween NuSTAR and Swift-XIIT is approximately 5% dMadsen et al. 

[2015J). 
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Table 3 

TVh values from different abundance models. 


Abund. a 

Model I 



Model II 



Model III 


Model 

AT H b x 2 c 

p -value 

7Vn b 

X 2c p -value 

A+ b 

x 2c 

p -value 

aneb 

1.06 ±0.02 2725.33 

5.5 x 

10“ 6 

0.41 ± 0.01 

2722.70 6.5 x 10-8 0.52 ± 0.01 2553.53 

1.9 x IO -2 

angr 

0.93 ± 0.02 2719.04 

8.1 x 

10-8 

0.37 ± 0.01 

2718.96 8.1 X 10 - 6 0.46 ±0.01 2551.61 

2.1 x 10“ 2 

f eld 

0.95 ± 0.02 2768.78 

3.7 x 

io - 7 

0.38 ±0.01 

2719.98 7.6 x 10-8 0.47 ±0.01 2567.20 

1.2 x 10“ 2 

grsa 

1.11 ±0.02 2753.07 

9.6 x 

io - 7 

0.44 ± 0.01 

2719.87 7.7 xl0 - 6 0.55 ± 0.01 2563.12 

1.4 x 10 “ 2 

wilm 

1.27 ±0.03 2826.47 

5.3 x 

10“ 9 

0.51 ± 0.01 

2722.94 6.4 xl0 - 6 0.64 ± 0.02 2588.34 

5.4 x 10 “ 3 

lodd 

1.32 ±0.03 2843.74 

1.4 x 

IO -9 

0.54 ± 0.01 

2718.83 8.2 x 10-8 0.67 ± 0.02 2592.26 

4.7 x IO -3 

aspl 

1.30 ±0.03 2827.83 

4.8 x 

IO -9 

0.52 ±0.01 

2722.14 6.7 x 10-8 0.65 ± 0.02 2584.44 

6.3 x IO -3 

a References 

— aneb: Anders & Ebihara 

1982) 

, angr: Anders & Grevesse (1989), feld: 

Feldman (19921), grsa: Grevesse 

|& Sauval (1998), wilm: Wilms et al. (^UUU 

, lod 

d: Lodders (2UU3f, aspl: jAsplund et al. (2 

009 



b In units of 10 22 cm 2 . 









c Each fit has 2408 degrees of freedom. 












Table 4 






Spectral fits to 0.5-79 keV phase-resolved Swift-X .RT and NuSTAR observations. 


Component Parameter 




Phase Range 






0 . 0 - 1.0 

0 . 0 - 0.2 

0.2—0.4 

0.4-0 .6 


0 . 6 - 0.8 

0 . 8 - 1.0 

const*tbabs*(bbody+powerlaw+powerlaw) 







const 

Cfpma 

0.981 ± 

0.015 

- 

- 

- 


- 

- 


Cfpmb 

0.977 ± 

0.015 

- 

- 

- 


- 

- 

tbabs 

JV H ( 10 22 cm -2 ) 

1.30 ± 

0.03 

- 

- 

- 


- 

- 

bbody 

k B T BB (keV) 

0.462 ± 

0.005 

0 474+ u - uus 

O /Iqq+0.009 

U.4oo_o 009 

0-467++ 9 

O /ic:o+0.008 

U.4oo_o Q08 

o 

“4 

OO 

1 + 

° P 

o o 

o o 

00 00 


norm a (IO -3 ) 

1.06 ± 

0.04 

0 qn + 006 
u - yu -0.06 

1 qk+ 0-07 
x - uo -0.06 

0.988+°® 

1 02 +0 ' 06 
l.UZ_o 06 

0 90+ b u 6 

u.yu—o.06 

powerlaw 

r s 

3.85 ± 

0.03 

q Of?+0.03 
^•^-0.03 

q Q.J+0.03 
d -° 4 -0.03 

q qi+0.03 
°- yi -0.03 

q Q7+0.03 
-0.03 

q Q9+0.03 
o,yz -0.03 


norm b 

0.18 ± 

0.01 

0 . 202 +;°°? 

0.192±g;°g| 

0 17q+0-°06 

u.i iy_ 0 _oo 6 

0 200~*~°' 007 
U.ZUU _ 0 Q07 

0 iq3 + 0 006 

u.iyo_ooo 6 

powerlaw 

r H 

0.29 ± 

0.05 

0 - 2 ±gd 

0-4±g;l 

0 3 + 0 - 1 


0 3 + 01 

0 3 + 01 
u - 0 -o.i 


norm b ( 10 “ 5 ) 

2.3 ± 

0.4 

o i + 0.8 

^-o.e 

4.6+8 

9 9 + 1-0 

z.z _ 0>7 


Z.D-o.7 

2 1+ 11 
z - 4 -0.7 

X 2 /dof 


2718.6/2408 

1184.8/1122 

1207.4/1120 

1011.0/1008 

1154.7/1045 

1070.2/1035 

p -value 


4.8 x 10 " 9 

9.4 x 10 ~ 2 

3.5 x IO ” 2 

4.7 x 10- 

1 

9.8 x IO -3 

2.2 x 10- 1 

const*tbabs*(nthcomp+powerlaw) 








const 

Cfpma 

1.001 ± 0.015 

- 

- 

- 


- 

- 


Cfpmb 

0.998 ±0.015 

- 

- 

- 


- 

- 

tbabs 

N b (10 22 cm -2 ) 

0.65 ± 0.02 

- 

- 

- 


- 

- 

nthcomp 

r s 

4.86 ± 0.04 

4 si +u.ut> 

4 75 +u,u7 

4 * '°-0.30 

4 89 +0 ' 08 

^•° y -0.15 

4 q 8 + b,us 
^• yo - 0.22 

4 94 4 " 9-10 
^•^— 0.21 


k B T BB (keV) 

0.346 ± 0.004 

0 344 + 0 003 
u.o^-0 009 

0 344 + 0 004 
u -°^-0.006 

0 34o + 0 004 
u.o-y:U_o Q04 

o 

CO 

CO 

OO 

1 + 
o o 

o o 
o o 
►p*. ^ 

0 346+ 0 003 

U.0^:0_o 004 


k B T e - (keV) 

> 37.3 

> 13.5 

> 4.8 

> 10.5 


> 6.7 

> 6.8 


norm b (IO -2 ) 

6.5 ± 0.2 

6.8 ± 0.1 

6.4 ±0.1 

6.0 ± 0.1 


6.8 ± 0.1 

6.4 ± 0.1 

powerlaw 

T/f 

0 75 + 0 - 05 
u -'°-0.04 

0.65++ 

0 79 “*"°' 1 
u -' y —0.08 

0 75+0 - 11 

0 71 + 01 
u - /J --0.09 

0-78l+ 9 


norm b ( 10 -4 ) 

1 1+0.2 
1 - 1 -o.i 

n q+0-3 
u - y - 0.2 

1 4+0-6 
^-O.S 

1 n+0.3 

J-^-o.s 


n q+O - 3 
u - y - 0.2 

i-o+l 

X 2 /dof 


2550.1/2408 

1190.4/1122 

1189.4/1120 

1012.0/1008 

1094.1/1045 

1058.8/1035 

p- value 


6.3 x 10~ 3 

7.6 x IO -2 

7.3 x IO" 2 

4.6 x 10- 

1 

1.4 x IO” 1 

3.0 x 10" 1 


a Normalization in units of L 39 /DJQ, where L 39 is the source luminosity in units of 10 39 ergs 1 and Dio is the distance to 
the source in units of 10 kpc. 
b Normalization in units of photons keV -1 cm - 


at 1 keV. 


instead of assuming that the j-bundle emerges from a po¬ 
lar cap, its footprint is allowed to have a ring shape. The 
assumption of axisymmetry reduces the number of free 
parameters and appears to be sufficient to fit the phase- 
resolved spectra. In future work, the energy-resolved 
pulse profiles could be included in the fit to constrain 
the axial distribution of the j-bundle. 

This more general model has the following parameters: 
(1) the power Lj of the e ± outflow along the j-bundle, (2) 
the angle a mag between the rotation axis and the mag¬ 
netic axis, (3) the angle /? 0 bs between the rotation axis 
and the observer’s line of sight, (4) the angular position 
9j of the j-bundle footprint, and (5) the angular width 
A9j of the j-bundle footprint. In addition, the reference 


point of the rotational phase, (f>o , is a free parameter, 
since we fit the phase-resolved spectra. 


We follow the method presented in Hascoet et al. 
(20141, and explore the whole parameter space by fit¬ 
ting the phase-averaged spectrum of the total emission 
(pulsed+unpulsed) and phase-resolved spectra of the 
pulsed emission. In order to get sufficient photon statis¬ 
tics, we used only three phase bins: ‘A’ (0.05-0.35), ‘B’ 
(0.35-0.70), and ‘C’ (0.70-1.05), roughly covering the 
primary pulse peak, the minima and the sub-peak, re¬ 
spectively. The bins are indicated in the last panel of 
Figure |T] The phase bin with the lowest flux is assumed 
to represent the “DC” (unpulsed) component; its spec¬ 
trum is subtracted from the total spectrum in each phase 
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Phase Phase 


Figure 7. Variation of nthcomp+PL model parameters as a function of rotational phase. In each plot, the solid black line shows the 
parameter value in each phase bin, dark and light gray regions show the 1-cr and 3 -<r error ranges respectively. The phase range is repeated 
twice for clarity. The dashed black line and dotted black lines show the value of the parameter in the phase-averaged spectral fit and the 
corresponding 3-cr error bars. Starting from the top left to bottom right, the plots show T#, power law normalization, FbTbbj r s an d 
nthcomp normalization, respectively. 


bin to obtain the spectrum of the pulsed component. The 
NuSTAR data are fitted above 16keV, where the hard 
component becomes dominant and the coronal outflow 
model has to account for most of the X-ray emission. 

The left panel of Figure [8] shows the map of p- values in 
the plane (a mag ,/3 0 bs)- The parameter space appears to 
be largel y degen e rate. For comparison with the res ults 


4.3 


uen 


of Hascoet et al. (2014| (discussed further in Section 
below), we also show the resulting p-value map w 
the footprint width is fixed to be A0j = 0 j/2, i.e. thin 
rings are excluded. Then the degeneracy of the param¬ 
eter space is significantly reduced, and the results are 
consistent with those of Hascoet et al. (2014). 

Using the obtained best-fit model for the hard X-ray 
component, we have investigated the remaining soft X- 
ray com p onen t. The procedure is similar to that in 
Hascoet et al. (20141: we freeze the best-fit parameters of 
the outflow model, and fit the spectrum in the 0.5-79 keV 
band including the Swift-X RT data. As in Hascoet et al. 
(2014), we find that the spectrum is well fitted by~~the 
sum of one blackbody, one modified blackbodjJ^j and 
the coronal outflow emission (which dominates above 
lOkeV). The (cold) blackbody and the (hot) modified 
blackbody have luminosities L c = 2.5(3) x 10 35 ergs -1 , 
Lh = 3.33(4) x 10 34 ergs _1 and temperatures kT c = 
0.408(3) keV, kT h = 0.85(1) keV similar to those fit by- 
model II in Section [375] The power-law tail of the modi¬ 
fied hot blackbody starts at E ta u = 5.7(1) keV. 

4. DISCUSSION AND CONCLUSION 


16 In this model, dubbed -B-Btail in|Vogel et al.|(:20li||, the Wien 
tail of the blackbody is replaced by a power-law “smoothly” con¬ 
nected at the photon energy -Ftaii- Here “smoothly” means that 
the photon spectrum and its derivative are continuous at Frail- 


We have described timing and spectral analysis of si¬ 
multaneous 0.3-79 keV Swift-XRT and NuSTAR obser¬ 
vations of 4U 0142+61. Using Fourier analysis we present 
the variation in pulse shape and pulse fraction over the 
soft X-ray and hard X-ray bands. We find a signif¬ 
icant change in pulse structure at the cross-over be¬ 
tween the soft-energy peak where the modified black¬ 
body emission is dominant and the hard-energy peak, 
where the magnetospheric tail emission is dominant. We 
do not find evidence for phase modulation in the 15 - 
40 keV lightcurve as reported by Makishima et al. (2014). 
We find that the phase-averaged spectrum is best mod- 
eled by a phenomenological nthcomp+PL model. The 
BB+2PL and 2BB+PL models that were traditionally 
used to fit the data do not provide statistically accept¬ 
able descriptions. Fitting the phase-resolved Swift-XRT 
and NuSTAR spectra of 4U 0142+61, we find that the 
spectral shape parameters do not show statistically sig¬ 
nificant variations compared to the phase-averaged fits. 
However, the normalizations of the spectral components 
vary significantly at phases corresponding to peaks and 
dips in the pulse profiles. Finally, we place constraints on 
the geometry of 4 U 0142+61 us i ng the electron-position 
outflow models of Beloborodov (2013a). 


4.1. Timing Analysis 

The low-energy pulse shapes measured from Swift- 
XRT and NuSTAR agree with the measurements of dH08 
gathered with XMM-Newton. In particular, the pulse 
profiles bear remarkable similarity with the data gath¬ 
ered on 2004 July 25 (dataset C) and on 2004 March 1 
(dataset B) and are less similar to the previous obser¬ 
vations (dataset A, gathered on 2003 January 04). The 
separation between the peaks in Swift-XRT (A(f = 0.35, 


















































































































p -value 


NuSTAR and Swift-XMT Observations of 4U 0142+61 


11 



Figure 8. Maps of p -values for the fit of the hard X-ray component with the coronal outflow model; the p -values are shown in the plane of 
(cumag, /5 0 bs) and maximized over the other parameters. The a ma g axis is common for all the plots. The p -value scale is shown on the left. 
The hatched green regions have p -values smaller than 0.001; the white regions have p -values greater than 0.1. Interchanging the values of 
a mag and /3 0 bs does not change the model spectrum, as long as the j-bundle is assumed to be axisymmetric. Therefore, the map of p-values 
is symmetric about the line of /5 0 t» s = a mag . Left: p -value map when A0j is thawed as a free parameter. Middle: p -value map when the 
footprint width is frozen to A0j = 0j/2. Right: p -value map when the footprint area of the ^‘-bundle, Aj, is restricted to be in the interval 
2.5 X 10 -3 < Aj/A^s < 10“ 2 (see discussion). 


in 0.3-1.5keV) matches that in the XMM-Newton data 
(A0 = 0.35, in 0.8-2.0keV). The separation between the 
dips (A(/> = 0.6) and the relative pulse heights also match 
well between the two data sets. Similarly, the pulse 
shapes and relative heights between NuSTAR 3-5 keV 
and 5-8 keV profiles and the XMM-Newton 2-8 keV pro¬ 
files are morphologically similar. 

Similarly, the NuSTAR 3-5 keV pr ofile agrees wit h the 
2-4keV RXTE pulse profiles of Dib et al. (2007) ob¬ 
tained between 2005 March and 2006 February. How¬ 
ever, there are increasing differences between the NuS¬ 
TAR profile and the RXTE pulse profiles at epochs go¬ 
ing backward from 2005 to 1996. The 6-8 keV profile 
between March 2005 and February 2006 shows a slightly 
broader main peak than the 5-8 keV NuSTAR pulse pro¬ 
files. Our low-energy pulse shapes agree well with 0.5- 
2keV and 2-10 keV XMM-Newton pulse profiles of |Gon-| 
zalez et al. (2010) obtained between 2006 July and 2008 
March with the"match improving as the compared epochs 
become closer. 

We find slight differences between the NuSTAR 20- 
35keV and 35-50 keV profiles and INTEGRAL 20- 
50keV profiles described in dH08. The NuSTAR profiles 
show a primary peak (at (f> = 0.2) that is 40% higher than 
the secondary peak (at </> = 0.7). In the INTEGRAL 
profiles, the peak separations are similar (A cj> = 0.53) 
but the peak count rates were equal. Similarly, we find 
that the NuSTAR 50-79 keV profiles show evidence of 
a double-peaked structure, with two sharp peaks sepa¬ 
rated by A <f> = 0.3. The corresponding 50-160 keV IN¬ 
TEGRAL profile shows a single-peaked structure. While 
it is possible that the pulse profile has changed, note 
that the 50-79 keV band would contribute only 35% of 
the photon flux as compared to the 50-160 keV energy 
band for a power-law spectrum with T = 0.65. Hence 
the difference in pulse profile may also be attributable 
to the difference in energy ranges. We also observe that 


the relative height of the primary pulse (at <f> ~ 0.3) 
compared to the pulse at <f> ~ 0.8 is decreasing with in¬ 
creasing energy through the 20-35 keV, 35-50 keV and 
50-79 keV plots. Hence it is not inconceivable that at 
energies higher than 79keV, the pulse at $ ~ 0.8 starts 
to dominate the pulse profile. 

4.1.1. Non-Detection of Precession 


After repeating the analysis steps of Makishima et al. 
(|20l4|) on 15-40 keV NuSTAR data of 4U 0142+61, we 
did not detect any phase modulation that can be in¬ 
terpreted as precession of the neutron star. The pulse 
profile from the NuSTAR data, while very consistent in 
shape and amplitude with the double-peaked profiles of 
dH08, are very different in shape and amplitude from 
the triple-peaked profiles of ME14 obtained after phase- 
demodulation. It is possible that the precession signal 
may be time-varying, having been detected in 2009 but 
not in 2007 and 2014. However, considering the necessary 
reconfiguration in the neutron star moments of inertia 
(A/ /I ~ 10” 4 ) and the corresponding reconfiguration of 
a 10 16 G toroidal magnetic field, it is surprising that the 
timing ephemeris, rotational spin-down and pulse pro¬ 
files remain consistent between 2007 and 2014. Further 
searches of phase modulation will help to confirm and 
understand the mechanics of this result. 

4.1.2. Comparisons with Other Magnetars 

The trend of pulse fraction as a function of energy 
varies from magnetar to magnetar though many have 
a p ulse fraction increa s ing with energy (see for ex am- 
ple Kuiper et al. 12006 den Hartog et al. p008 bfla| ). In 
dUOTl^TST^eoDserve tHaO T FRMFlncreases up to a 
value of 20% and possibly shows a small decline to¬ 
wards 40keV or possibly stays constant at «20%. In 
IE 2259+596, PPrms was seen t° monotonically rise to 
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approximately 70% at 20keV (Vogel et al.||2014 | and in 
IE 1841—045, PErms was seen to rise to a value of ap¬ 
proximately 17% at lOkeV, decrease to 12% at 20keV 
and ris e again to approxim ately 17-20% between 30- 
79keV ([An et al.|2013||2015|). At the same time, PF area 
was measured to increase from 25% at 1-2 keV and in¬ 


crease to 50% at 50 keV. For 1RXS J170849-400910 Aden 


Hartog et al. 2008a), the pulse fraction (reported as 
PT'area, not PT'rms] was shown to be nearly constant at 
approximately 40% between an energy range from 0.7 
200 keV. In fact, the pulse fraction decreases slightly from 
about 50% at 1 keV to about 30% at 3 keV, rising back 
to about 40% at higher energies. 

The X-ray pulse profiles of magnetars, affected by the 
geometry of the magnetic field and rotation axis, are sim¬ 
ilarly diverse. The pulse profiles of 4U 0142+61 are pri¬ 
marily double peaked for most energy bands with each 
pulse width being 60 « 0.25. Compared to these, the 
pulse profile of IE 1841—045 (An et al. 2013. 2015) is 
comprised of large, single-peaked - humps that are about 
60 « 0.75 wide (except for the double-peaked struc¬ 
ture emerging between 23.8 -35.2 keV). The pulse p rofiles 


of 1R.XS J170849—400910 (den Hartog et al. 


2008a) are 


dominated by a single pulse peak with a width d0 0.35 
at most energy bands; however, there is a distinct shift 
between pulse positions below and above 8 keV, suggest¬ 
ing that the dominant emission mechanism changes dras¬ 
tically. The pulse profiles of IE 2259+586 show compli¬ 
cated structure, with narrow peaks ( 60 ss 0 .25) that can 
possibly shift slightly with energy (Vogel et al. 2014). 
The pulse profiles of 4U 0142+61 are therefore morpho¬ 
logically more similar to those of 1RXS J170849—400910 
than those from IE 2259+586 and IE 1841—045. 

A possible source for these differences may be the 
size and geometry of the hot-spot emitting area on 
each magnetar. Comparing the size of the j-bundle, 
9j in the outflow model fits (see Section 4.3) suggests 
a rough pattern, albeit in a very limited sample size. 
For IE 1841—045, the magnetar with the broa dest pulse 
profiles, 0j < 0.4rad (An et al. 2013 2015) while for 
1RXS J170849—400910 and 4U 0142+61 with narrower 
pulse profil es, 61,- < 0. 15 r ad and < 0.23 rad, respec¬ 
tively (Hascoet et al. 2014). The outflow model fit for 
IE 2259+586, which shows a complicated narrow pulse 
profile, statistically prefers a complicated ring-shaped j- 
bundle with 0.4 rad < 9j < 0. 75 rad an d A 9j/6j, the 
ring-width fracion, < 0.2 (Vogel et al.||2014|. 


4.2. Spectral Analysis 

We fit different spectral models to the soft and hard 
energy spectra from Swift-XRT and NuSTAR. We find 
that the 2BB+PL and BB+2PL models do not fit the 
cross-over region (approximately 5-15 keV) of the spec¬ 
trum well. This causes a distortion in the fitting of the 
hard power-law and a residual is left at the high ener¬ 
gies (> 50keV). The nthcomp+PL model provides a sta¬ 
tistically better fit than the two models, especially for 
the cross-over region. The hard PL index Th measured 
from this model best matches the Th = 0.65 + 0.09 mea¬ 
sured after restricting the energy range to be between 
20-79 keV. 

We find that the spectral turnover Ps — Th = 3.56 
(BB+2PL model) and Tg — Th = 4.11 (nthcomp+PL) 
are higher than the values reported for the total flux 


Ts — = 2.6 reported by Kaspi & Boydstun (2010). 

Using the independent value of L'h = 0.65 + 0.09 slightly 
increases the discrepancy. Placing these v alues on the 
r, 9 — Th vs log(B/10 14 G) plot (as shown in Vogel et al. 
20141, does not change the observed decreasing trend 
between Pg — T# and log(R/10 14 G). 

Fitting the same models to the phase-resolved spectra 
shows that the spectral shape parameters ( Th , A+Tbb 
and r s) are consistent within 3-cr error bars to the val¬ 
ues measured from phase-averaged spectral fits. How¬ 
ever, the normalization of the hard X-ray and soft X-ray 
components varies significantly as a function of phase. 
From Figure [7] we can identify the increase in the hard 
power-law normalization in the 0.2-0.4 phase range with 
the peak in the 20-35 keV and 35-50 keV pulse profiles 
and the dip in the soft X-ray components normalization 
(nthcomp or soft power law, depending on the model fits) 
with the dip at phase 0 = 0.5 in the 3-5 keV and 5-8 keV 
bands. This suggests a clear differentiation between the 
low-energy and high-energy spectral components. 

We find that the best-fitting nthcomp+PL model yields 
IVh = (6.5 ± 0.2) x 10 21 cm' 2 , consistent with that mea¬ 
sured by D06b and also co nsistent with later broadband 
fits by Enoto et al. (2011). 


4.3. Outflow Model 

We find that the coronal outflow model provides con¬ 


4U 0142+61. 

Hascoet et al. 

(2014 

) obtained a sirni- 

lar conclusion 

by fitting the lNl'EGl 

RAL phase-resolved 

spectra of 

den Hartog et al. 

(2008b 

). Their model as- 


emerging from a polar cap on the star or a thick ring 
A0j = 0j/2. An excellent fit was provided by this model 
in a small region of parameter space, giving strong con¬ 
straints on a mag and / 3 0 hs- Motivat ed by the recent anal¬ 


ysis of 1E2259+586 (Vogel et al. 2014), we explored a 
more general outflow model that allows the footprint of 
j-bundle to be a ring of arbitrary thickness A0j. We 
found that a thin-ring configuration is also able to fit 
the phase-resolved spectrum of 4U 0142+61, and in a 
broader range of parameters. This degeneracy is absent 
in IE 1841 —045 , where on ly a thick ring or a polar cap 
is allowed (jAn et al. 2015). 

The soft X-ray component (below ~ 10 keV) is well fit¬ 
ted by the sum of one cold blackbody and one modified 
hot blackbody. The cold blackbody covers a large frac¬ 
tion of the neutron star area, A c ~ 0.7.4ns- The emis¬ 
sion area of the hot blackbody is small, Ah ~ 0.005 Ans- 
In the coronal outflow model, the footprint of the j- 
bundle is expected to form a hot spot, as some particles 
accelerated in the j-bundle flow back to the neutron star 
and bombard its surface. If the modified hot blackbody is 
interpreted as the thermal emission from the footprint, 
then the measured Ah can be used as a constraint on 
the footprint of the active j-bundle. The right panel of 
Figure [5j shows the NuSTAR p-value when the j-bundle 
footprint area, Aj = n sin 2 9 } , is restricted to be between 
Ah/2 = 2.5 x 10 3 x4ns and Ah x 2 = 10 2 Ans- Then 
the degeneracy of the model is reduced and a broad re¬ 
gion of the parameter space around the line a mag = /3 0 bs 
becomes excluded. 

The outflow model predicts the X-ray flux below 
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~1 MeV to be dominated by photons polarized perpen¬ 
dicular to the magnetic field while an excess of parallel- 
polarized photons i s expected th rough photon splitting 
at higher energies (Beloborodov 2013aI. The model of 
magnetospheric emission will be crucially tested by fu¬ 
ture X-ray polarime try instruments such as ASTRO-H- 
SGD (Uyimaetal. 20101 during flares), ASTR OSAT- 


CZTI (Uhattopadhyay et al. 2014|), POLAR (Produit 
et al.(20L)5[ > and X-Calibur ( jBeilicke et al.|2014 |. 


In conclusion, we have presented a timing and spec¬ 
tral analysis of simultaneous 0.5-79 keV observations of 
4U 0142+61 using Swift-X RT and NuSTAR. The rota¬ 
tional period of 4U 0142+61 is consistent with that ex¬ 
pected from extrapola tion of t he timing solution since the 
last glitch (Dib & Kaspi 2014). We have not detected the 
55-ks time period, 0.7-s amplitude phase modulation in 
the 15-40 keV Suzaku- HXD data from 2007 reported by 
Makishima et al. (2014) that was ascribed to the free- 
precession of 4U 0142+61. While this precession may be 
time-varying, the consistency of the rotational ephemeris 
and pulse profile between 2007 and 2014 needs to be ex¬ 
plained. We have shown that the pulse profile changes 
character (dominance of the first harmonic vs the second 
harmonic) at around 30 keV. While the low-energy pulse 
profiles were consistent with pre viously prese nted pulse 
profiles (b etween 2006 and 2008: Dib et al.|2007 |Gonza- 
lez et al. 2010), we have observed morphological differ¬ 
ences between the hard energy pulse profiles of NuSTAR 
and INTEGRAL. We have shown that the RMS pulse 
fraction has an increasing trend with energy, reaching 
a value of up to 20%, however, it shows some evidence 
of a decrease at about 40keV similar to that observed 
in IE 1841—045 and contrary to the smooth increase of 
pulse fraction in IE 2259+586 that increases to nearly 
80%. 

We have shown that the energy spectrum of 
4U 0142+61 between 0.5-79 keV is better described by a 
Comptonized blackbody + hard PL model than the pre¬ 
viously used BB+2PL or 2BB+PL models with a hard 
power-law (T# = 0.65 ± 0.09) dominating the spectrum 
above 20 keV. The low-energy spectrum (< 10 keV) may 
still be fit with the BB+PL model, however, this model 
cannot fit the observed spectrum between 10-20 keV. 

We have fitted the phase-resolved spectra of 
4U 0142+61 with the e ± outflow model of Beloborodov 


( 201 3b) using the analysis method of Hascoet et al. 
(2014|). Our results show that the outflow model gives 
a consistent physical description of the phase-resolved 
spectra, and the results are consistent with those derived 
from INTEGRAL data. We found that significant degen¬ 
eracy appears in the inferred parameters of the inclined 
rotator a mag and /3 0 bs if the footprint of the j-bundle is 
allowed to be a thin ring. The degeneracy is significantly 
reduced if the footprint area A, is restricted to be similar 
to the area of the blackbody hotspot that covers 0.5% of 
the neutron star surface. 
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